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Abstract. We review chiral (Klein) tunneling in single-layer and bilayer graphene 
and present its semiclassical theory, including the Berry phase and the Maslov index. 
Peculiarities of the chiral tunneling are naturally explained in terms of classical 
phase space. In a one-dimensional geometry we reduced the original Dirac equation, 
describing the dynamics of charge carriers in the single layer graphene, to an effective 
Schrodinger equation with a complex potential. This allowed us to study tunneling 
in details and obtain analytic formulas. Our predictions are compared with numerical 
results. We have also demonstrated that, for the case of asymmetric n-p-n junction 
in single layer graphene, there is total transmission for normal incidence only, side 
resonances are suppressed. 



1. Introduction 

Since this paper is prepared for the proceedings of the Nobel symposium on graphene we 
do not start with general explanations what graphene is and why it is important, it will 
be very well described in other presentations. We just refer to reviews [H |2l [3l HI |5l [6l [7] . 
Our particular subject is chiral, or Klein (as it was called in [8]) tunneling. This is one 
of the key phenomena determining peculiar electronic properties of graphene. In light of 
possible applications, the Klein tunneling protects high charge carrier mobility despite 
unavoidable inhomogeneities. At the same time, due to the Klein tunneling graphene 
electronics cannot copy the standard semiconductor one: if you make graphene transistor 
based on n-p-n junction just like for silicon, it will not be efficient since you will not be 
able to lock it. These two remarks illustrate the importance of the subject, the more 
detailed discussion is presented below. 

The paper consists of two pieces. The first one (Sections [2] - [5]) preserves the 
historical line of thoughts and presents the motivation of the problem, from [9] to [8]. 
In the second part (Sections [6] -[12]) we present a systematic semiclassical theory of the 
chiral tunneling, together with numerical results. 

2. The Klein paradox 

Soon after the discovery of the Dirac equation, O. Klein \^ noticed one of its strange 
properties which was afterwards called the "Klein paradox". Klein considered the 
original four by four Dirac equation, which governs the dynamics of a spin one half 
particle moving in three-dimensional space. To make a direct connection to the case of 
graphene without changing the essence of the paradox, we will consider a two by two 
matrix equation for a particle propagating in two-dimensional space: 



Here m is the mass of the particle, c is the speed of light and u{x,y) is the potential 
energy. 



where = {ipi,ip2) and the Hamiltonian 



(1) 




(2) 
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To demonstrate the essence of the paradox we consider normal incidence on a 
one-dimensional potential barrier, which means that u = u{x) and tpi = ipi{x). Then 
equation ([T]) takes the form 

dijj2 



^he^ = {E- me- - u{x)) , 



in,e— — = iE + me — u{x)] il)2 ■ 



dx 

To make the problem exactly solvable we use a step-wise potential 



j 0, X < 0, 

u[x) = { ^ (4) 

where Uq is a positive constant. We consider a general scattering problem with an 
incoming wave '^in{x) and a reflected wave '^r{x) for x < 0, 

^(x) = *i„(x) +r*^(x) , (5) 

and a transmitted wave \l/t(x) for x > 0, 

$(x) = t^t(x) . (6) 

The x-dependence of the solutions for x < can be written as exp(±iA;x), where 
the wave vector k satisfies the relativistic dispersion relation E'^ = fi^e^k'^ + m^c^ as can 
be found by diagonalizing equation ([3]) with m = 0. Alternatively the wave vector can 
be written as 



m^e 



2^4 



One easily sees that there are three distinct regimes, two of which are classically allowed, 
namely E > mc^ corresponding to electron states and E < —me^ corresponding to hole 
or positron states. There is also a classically forbidden region —me^ < E < me^ where 
the wave vector k is imaginary and we have evanescent waves. In what follows we will 
assume that we are in the electron regime. By calculating eigenvectors of equation ([3]) 
one obtains for the wavefunctions to the left of the barrier 

'^^n{x) = ( M e*'" (8) 



and 



where 



vl/,(x) = f M e-'^^ , (9) 



—a 



E — mc^ , , 

a = \- -. 10 

\ E + me^ ^ ' 

To the right of the barrier we have a new wave vector g, which satisfies the 
relativistic dispersion relation {E — uqY = h'^c-q'^ + m^c^, or 



'(uo - EY - m^e^ 
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Consider a jump 

Uo>E + mc^, (12) 

since in this case the paradox arises. The wave vector q is real and we have a propagating 
wave on the right side of the barrier. Note however that this particle belongs to 
the hole continuum rather than to the electron one. For smaller values of Uq, there 
are either propagating electrons on both the left and the right side of the barrier, 
when uq < E — mc^, or evanescent waves on the right side of the barrier, when 
E — mc^ < Uq < E + mc^. Solving the Dirac equation ([3]) on the right side of the 
barrier, one obtains for the transmitted wave 




where 



= si ° z.:z, - (14) 



Iuq — E — mc^ 
Uq — E + mc^ 
From the continuity of the wavefunction at x = 0, 

^in + r^r\x=-o = ^tU=+o , (15) 



we find 



a/3 + 1 



(16) 



a/3-1 ■ 

For the considered case we have < a, /3 < 1, so that r < and 

''HH'^^(i±^y>i. (1.) 

To treat reflection and transmission coefficients properly one has to look at the 
probability current density for the one-dimensional Dirac equation 

= c¥a^^ = c{i)li)2 + ^2>i) , (18) 

which is a conserved quantity. When we look at the current density (fTSj) we see that it 
takes values 2ac for the incoming wave and — 2aci? for the reflected wave. Therefore R is 
nothing but the reflection coefficient and we come to the conclusion that the amplitude 
of the reflected wave is larger than the amplitude of the incident one. This strange effect 
that occurs when condition ([T^ is fulfllled was initially called the Klein paradox. In our 
further discussion we will follow [10] and [11]. For a rather complete list of references 
see [I2]. 

First of all note that the current density (|T8l) on the right hand side equals — 2|tp//3, 
indicating that there is something wrong with the deflnition of the transmitted wave. 
What exactly is wrong was pointed out by Pauli, who noticed that the group velocity 
for the case of equation (fT2l) . 
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is opposite to the direction of the wave vector q. Since the group velocity determines 
the direction of propagation, the transmitted wave f|T3l) corresponds (for positive q) to 
a particle moving to the left instead of to the right. Therefore we should define our 
outgoing wave as 

^t{x) = ( ^^^^ ) e-"- , (20) 
which gives the currenty density 2|tp//3. When we once again calculate r from 



equation ( 1T51) . it is seen that 

R=\r\^=(\^y <l. (21) 
^ ^ \l + a/3j ^ ^ 

which is always smaller than one. Therefore the formal paradox disappears, see also [T3] . 

The paradox reappears when we consider the problem from a different angle. 
Instead of an infinitely broad barrier we will consider a finite barrier, 

"W^r I:!:: 

The problem with the choice of the transmitted wave on the right side of the barrier has 
now disappeared, since it is simply t'^in- Within the barrier one now has to consider both 
modes exp(±iga;), representing the most general solution. Reflection and transmission 
coefficients are then obtained from the continuity of the wave function at x = —a and 
X = a, which gives after straightforward calculations (see e.g. [10] and [H]) 

^ il-a'/3rsm'{2qa) 

4a^/3^ + (1 - a2/32)2 sm\2qa) ' ^ ' 



T = ^ (24) 



4a2/32 + (i_c,2/32)2sin^(2ga)' 

There is no paradox in these expressions, since 0<i?<l,0<T<l and i? + T = 1 as 

it should be. Note that we have total transmission through the barrier when 

iVvr , , 

ga = — , (25) 

with integer A^. 

We can consider an infinitely broad barrier by letting a go to infinity in the above 
expressions. As a becomes very large while other parameters remain fixed the sine 
will oscillate very rapidly. We can then average over the fast oscillations and replace 
sin^(2ga) by its average value ^ to obtain the expressions 

= 8a2^2 + (1 _ c^2^2)2 (26) 
8a2/32 (27) 



^ 8a2/32 + (1 _ ct2^2)2 

One may be surprised that the results ( 1211) and fl2B]) do not coincide. It is however well 
known from electromagnetic wave theory [15] that the refiection coefficients for the two 
situations should differ. 
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From the last result we see once again that the paradox has disappeared in its 
mathematical form, but has reappeared as physically counterintuitive behaviour. In 
non-relativistic quantum mechanics a particle can tunnel through a classically forbidden 
region E < u{x), but the probability is exponentially small when the barrier is high and 
broad. In the semiclassical approximation the transmission through the barrier with 
turning points 0:1^25 which satisfy E = u{xi^2), is given by 

2 /"^a 



T = exp (^-- y dx,j2m{u{x) - E)j , (28) 

where m is the mass of the particle. For a relativistic particle incident on a sufficiently 
high barrier ( !T2|) the situation is dramatically different. In the limit a — t- 00 the 
probability of penetration (127|) is in general not small at all. Even for an infinitely 
high barrier {uq — ?■ 00) one has /3 = 1 and 

E'^ — w?c^ 
E^ — ^m^cr 

This is of the order of one when E — mc^ is of the order of mc^, while it is approximately 
equal to one in the ultrarelativistic limit 

E > mc^ . (30) 

This is the contemporary formulation of the Klein paradox [lOj; quantum relativistic 
particles can tunnel with large enough probabilities through barriers of arbitrarily large 
height and width. 

The tunneling effect can be hand-wavingly explained with the help of the Heisenberg 
uncertainty principle. Since one cannot know both momentum and position with 
an arbitrary accuracy at a given instant, one cannot separate the total energy into 
a potential and a kinetic part. So the kinetic energy can be "a bit" negative. In 
the relativistic regime the restriction is much stronger [16]: one cannot even know 
the coordinate with an accuracy higher than hc/E. Therefore relativistic quantum 
mechanics cannot be mechanics, but can only be field theory [17]. This theory will 
always contain particles and antiparticles and to measure the coordinate better than 
hc/E one needs to apply such a high energy that particle-antiparticle pairs will be 
created. The original particle whose coordinate one wanted to measure will then be 
lost among the newly-born particles. A full field theoretic treatment of the problem 
was given in Ref. [18]. The most important point is that although the problem of a 
high enough barrier looks like a static problem, this is actually not the case. One needs 
to study carefully how the state is reached and this involves positron emission by the 
growing barrier. For a more detailed discussion of the role of electron-positron pairs in 
the Klein paradox, see [T9] . 
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3. Klein tunneling in single layer graphene 

The Hamiltonian for charge carriers in graphene near conical points K and K' is given 
by the massless Dirac Hamiltonian 

H = V (a^p^ + ayPy) + u{x, y) , (31) 

where V is the Fermi velocity V ~ c/300. To consider normal incidence on the one- 
dimensional potential barrier ( 122|) in this case, we can simply put m = in our previous 
results. From equations ( JTOl) and it is seen that a = [5 = 1. Therefore we have 
T = 1 and i? = in equations fl26l) and fl27|) . regardless of the height of the potential. 
This result is not related to the specific shape of the potential |20] . 

This property has an analog in two and three dimensions with u = u{x,y) or 
u = u{x,y,z), namely that backscattering is forbidden. This was found long ago 
for scattering of ultrarelativistic particles in three dimensions (see [211 IE])- An 
important consequence of this property for carbon materials was noticed in |20j . 
Absence of backscattering explains the existence of conducting channels in metallic 
carbon nanotubes, while in a non-relativistic one-dimensional system an arbitrarily small 
disorder leads to localization [22j . 

The consideration in [20] is very instructive since it explicitly shows the role of 
the Berry phase and time-reversal symmetry, but it is also quite cumbersome. Here 
we present a somewhat simplified scheme of this proof. To this aim we consider the 
equation for the T-matrix (see e.g. [23] ) 

f = u + uGof , (32) 

where u is the operator corresponding to the scattering potential. 

Go = lim 1 , (33) 

E-Ho + i5 

is the Green's function of the unperturbed Hamiltonian Hq and E is the electron energy, 
which is assumed to be larger than zero. If Hq is the Dirac Hamiltonian for massless 
Dirac fermions (13T|) . we have 

Goir, r') = I ^Go(q) exp[^q(r - r')] , (34) 

where 

A / ^ _ 1 _ 1 e + qa- 

^0^^^ - E- nVqa- + iS~ hV{e + i5f - ' ^^^^ 
with e = E/hV. The probability of backscattering can be found by iterating 
equation ( 132|) and is proportional to 

T(-k, k) = (-k \u + uGqU + uGouGqU + . . . I k) = T(^) + r(2) + . . . , (36) 
where T*^") is the contribution proportional to u^. 
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We can always choose axes such that k || Ox. In this case |k) and |— k) have spinor 
structures ^ | ^ \ ) respectively. Therefore, if T is the two by two matrix 

f = To + To- , (37) 

one has 

T(-k, k) ~ T,(-k, k) + iTy{-k, k) (38) 

Now keeping in mind that V is proportional to the identity matrix one can prove term 
by term that all contributions to Tj^(— k, k) and T;j(— k, k) vanish by symmetry. Actually 
this is because T(k) ~ k || Ox; from the vectors k and — k one cannot construct anything 
with nonzero y or z components. Strictly speaking this argument is only enough for an 
isotropic potential; for a generic case one has to do a term by term analysis based on 
expansion ( |36l) . see Ref. [20]. For two nonparallel vectors ki and k2 one can construct 
a matrix with nonzero y or z components, since one of the vectors has a nonzero y 
component, so that ki x k2 || Oz. 

When one thinks about electrons in quantum electrodynamics, it is not easy to 
create potential jumps larger than 2mc^ ~ 1 MeV. Similar phenomena take place in 
electric or gravitational fields ([2ll ES]; see [I2] for a detailed list of references), but 
the context is always quite exotic, such as collisions of ultraheavy ions or even black 
hole evaporation. There were no experimental data available which would require the 
Klein paradox for their explanation. However shortly after the discovery of graphene 
it was realized that Klein tunneling is one of the crucial phenomena for graphene 
physics and electronics ^ . Soon after this theoretial prediction the effect was confirmed 
experimentally [26], [27] . 

Considering possible applications, Klein tunneling in graphene is rather bad news. If 
one copied the construction from a silicon transistor to graphene, it would be impossible 
to lock the transistor. One would need to open a gap in the spectrum to be able to lock 
it. At the same time it is good news as well: due to the Klein paradox inhomogeneities in 
the electron density do not lead to localization and their effect on the electron mobility 
is not very essential [8]. 



4. Tunneling trough a stepwise barrier 

Let us now consider a massless Dirac fermion incident on the potential barrier fl22|) with 
positive energy under an angle (p, as it was done first in [8]. Of course, the potential 
cannot be sharp on the atomic scale, since this would induce Umklapp scattering 
between different valleys. Therefore by a step-wise potential we mean that the electron 
wavelength is much larger than the typical spatial scale of the potential /, which is 
in turn much larger than the size of the unit cell. 

Within this assumption the solution is each region is given by traveling waves 
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proportional to exp{±ikxx) exp{±ikyy) , where and ky satisfy the dispersion relation 

'-^f^k' = ikl + kl), (39) 

as can be found from equation (l3T|l . Similarly to the original Dirac equation we can 
distinguish three distinct regimes from this equation. For uq < E — hV\ky\ we have 
electrons and for uq > E + %V\ky\ we have holes, while the region E — hV\ky\ < uq < 
E + is classically forbidden. As was done for the case of the massive Dirac 

equation, we will now require that the potential Uq in equation fl22|) satisfies 

uo>E + hV\ky\, (40) 

so that we have hole states within the barrier. 

Let us denote by k the wave vector for |x| > a and by q the wave vector for |x| < a. 
At the potential jump the momentum in the y direction should be conserved, so that 
the new angle 6 is related to the new wave vector q by 

k sin (f) = ky = qy = q sin 6 . (41) 

From equation fl3Tl) we see that the second component of the wavefunction is related to 
the first by 

^P2 = sgn{E - u^)e^iP^ , (42) 
SO the solutions in the three regions are given by 

j ^ikxX ^ikyy _j_ ^ | j ^—ikxX ^ikyy ^ ^ ^ 



^{x,y) 



A I \g ) e''^-''e'''yy + B I ] ) e-'^-'' e'^yy , -a < x < a (43) 



s e \ —s e 



where we have introduced s = sgn{E), s' = sgn{E — uq), k^ = kcoscf) and q^ = qcos9. 
Note that the reflected particle moves under the angle vr — 0, assuming that the angle 
changes from —tt/2 to 37r/2, so that we have the phase — exp(— i0) for the reflected 
wave. We can now determine the reflection coefficient r, the transmission coefficient t 
and the coefficients A and B as before, from the requirement that the wavefunction is 
continuous at x = ±a. 

Finally the result is given by 

^ is-2ik n ■ N sin - ss' sin 6' 

r = 2e"^ sm(2g^a) — — — . . (44) 

^ ' ss' [e-^'i-" cos(0 + 9) + e^"?-"^ cos(0 - 9)] - 2i sm(2g^a) ^ ' 

For the case under consideration we have ss' = —1, since the signs of E and E — uq are 
opposite. The transmission probability can now easily be calculated as 

T = |tp = 1 - |r|2 . (45) 

From equation (H4l) we immediately see that the reflection is zero for normal incidence, as 
we proved for a more general potential in the previous section. There are also additional 
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angles, called "magic angles" , at which the reflection coefficient is zero and we have full 
transmission. They are given by the condition 

qxa = iv| , (46) 

where N is an integer. 

We can compare the behaviour of electrons in single layer graphene with the 
behaviour of normal electrons. When the potential barrier contains no electronic states, 
the transmission decays exponentially with increasing barrier width and height, [28J, so 
that the barrier would reflect electrons completely. But since single layer graphene is 
gapless, it seems more appropriate to compare it to a gapless semiconductor with non- 
chiral charge carriers, a situation which can be realized in certain heterostructures [29| 
130] . For this case we find 

^ 4A;^g^exp(2ig^a) 

{q + k^y exp{-2iq^a) - {q^ - k^y exp{2iq^a) ' 
where k^^ and q^ are the x-components of the wave vector outside and inside the barrier, 
respectively. As in the case of single layer graphene there are resonance conditions at 
which the barrier is transparent, given by 2qxa = Ntt, where is an integer. For 
normal incidence we see that the transmission coefficient is an oscillating function of 
the tunneling parameters and can exhibit any value between zero and one. This is in 
contrast to single layer graphene, where the transmission is always perfect. 

5. Klein tunneling in bilayer graphene 

Bilayer graphene consists of two layers of graphene on top of each other, the second 
layer being rotated by 120 degrees with respect to the ffist one. In this configuration 
the sublattices A he exactly on top of each other and the hopping parameter 71 between 
them is approximately 0.4 eV [311 [32] . while the in-plane hopping parameter 70 = t 
is approximately an order of magnitude larger. When we consider only low energy 
excitations, [E'l, |-E — mo| <^ 2|7i|, the effective Hamiltonian is given by [331 



H=( ° , - \ 

where the effective mass m = '-fi/2V^ ~ 0.054me, me being the free electron mass [35j . 
There is also hopping between the B sublattices of both layers, which is denoted by 
73 ~ 0.3 eV. When we include this parameter into the description an extra term is 
added to the Hamiltonian, which corresponds to so-called trigonal warping. This effect 
is however only important for small wave vectors ^35j , we will exclude it assuming that 
ka,qa > 737i/7o- 

Let us consider an electron incident on the potential step fl22l) under an angle 0, 
as was done in [8]. Since the potential is constant in the ^/-direction we can write the 
solution as 

<if{x,y) = ^{x)e'''yy . (49) 
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Inserting this into equation ([T]) with the Hamiltonian f l48p . we obtain 

{^-^..^(^^^'..^^',.. (50) 

The solutions are therefore given by propagating waves exp[±ikxx) and exponentially 
growing and decaying modes exp(± 
,2 ,2 2m\E — u\ 

kl + kl = , (51) 

2 2 2m\E — u\ 

f^x ~ = • (52) 

The presence of evanescent modes is markedly different from both the Schrodinger case 
and the Dirac case. Once again there are three regimes. There are electron states for 
Mo < E — h'^ky/(2m) and hole states for uq > E + h'^ky/(2m), while the region in between 
is classically forbidden. In what follows we assume that uq in equation ( l22l) satisfies 



uo>E + —^. (53) 
2m 

To find the spinors that are the solutions to equation (!50!) we note that the components 
are related by 

( d \\ 2m{E-u) , 

as can be seen from the Hamiltonian (1481) . 

Now let k = \/2mE/h be the wave vector for the propagating modes in the region 
|x| > a, while q = ^2m{uo — E)/h is the wave vector in the region < a. Then the 
solution for x < —a is given by 

= ( ) e''- + 6. ( ) + ( ) . (55) 



where ky = /csin0, kx = /ccos0, s = sgn(£'), = Jk'^ + 2ky = kJl + sin^ (f) and 



finally /ii = (y 1 + sin^ — sin0)^. The amplitude Oi is the amplitude for the incoming 
wave in this expression, while bi corresponds to the reflected wave. For a; > a we have 
the general solution 

*(-) = '.3( J.,)e*- + <,(_;^^Je-«, (56) 

where as is the transmission coefficient. Inside the barrier we need the most general 
solution with two propagating modes and two modes with real exponentials, 

where qy = qsinO = ky because the transverse momentum is conserved. Furthermore 
qx = qcosd, s' = sgn{E — uq), Xx = qVl + sin^ 6 and /12 = {Vl + sin^ 6 — sin6')^. 
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Now the coefficients Oj, bi, Ci and di have to be found from the continuity of ipi{x) 
and the derivative dipi/dx at the points x = ±a. When the problem is solved numerically, 
one sees that the transmission probability at normal incidence is exponentially small. 
Similar to the case of single layer graphene, there are once again "magic angles" in 
the spectrum, at which there is total transmission. The existence of magic angles in 
bilayer graphene has the same consequences as in single layer graphene, meaning that 
we cannot lock a conventional transistor made from bilayer graphene. 

For the case of normal incidence (p = 6 = we can also solve the problem 
analytically. The transmission coefficient is given by 
^ Aikq exp{2ika) 

{q + iky exp(— 2ga) — (g — ik^ exp(2ga) 
which is indeed exponentially small. When we let a go to infinity, the transmission 
probability T = |tp becomes zero at normal incidence. Furthermore for a single n-p 
junction with uq ^ E the following analytical solution can be found for any 

r=— sin2(20), (59) 

Mo 

which also gives T = at normal incidence, in contrast to the case of single layer 
graphene, where normally incident electrons are always transmitted. It is also different 
from the case of normal electrons, where the transmission is given by equation (1471) . 



6. Dimensionless variables and parameters 



In sections [HI [5] it was discussed that the wavefunctions \l/ of charge carriers in single 
layer and bilayer graphene in a one- dimensional geometry obey equations 





V 



Px - iPy 
Px + iPy 



+ u{x/l) - E 



^ = 0, 



and 







{Px - iPy) 



+ u{x/l) -E 



^ = 0, 



(60) 



(61) 



_2m \ (p^ + ipy^ 

respectively. Here / is a characteristic scale of a potential change. In dimensionless 
variables f l60l) takes the form 





+ u(x] 



E 



^ = 0, 



(62) 



Px - Wy 
Px + iPy 

where x = x/l, p^ = —ihd/dx, Py = Py/po, h = h/pol, u = u/Vpo and E = E/Vpo. We 
denote some characteristic value of \u — E\ as Vpo. 



Analogously, ( ]6T1) can be rewritten as 

{px - iPyf ^ 



{Px + iPy 







+ U[Xj 



E 



^ = 0, 



(63) 



with X = x/l, Px = —ihd/dx, Py = Py/po, h = h/p^l, u = 2mu/p1 and E = 2mE/p1. 
We denote some characteristic value of |m — E\ as p'^/2m. 
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H 



Thus we can introduce dimensionless Hamiltonians (we omitted tildes): 

\Px + iPy ; 

for a single layer and 



(64) 



(65) 



, iPx + iPyY 

for a bilayer. In both cases there are two substantial parameters in the problem: h and 

Py 



7. Standard semiclassical treatment 



Charge carriers in single layer graphene are described by the Hamiltonian flMl) . This 
Hamiltonian describes simultaneously coupled electron and hole states. According to 

can be diagonalized up to any order of 



Appendix A , in adiabatic approximation ([& 
h 1. The obtained scalar Hamiltonians describe electrons and holes separately. 
The diagonalization is based on a series of unitary transformations of the original 
Hamiltonian and traces back to the ideas of the Foldy-Wouthuysen transformation [36] 
and the Peierls substitution in Blount's treatment [3^. We use its variant [381 139] . 

Effective electron and hole Hamiltonians and L~ can be written as series with 
respect to the small parameter h\ 

L^{p^,x,h) = L^{p^,x) + hLf{p^,x) + h^L^{p^,x) + ... (66) 

To be precise we will assume that any function of px and x is defined in such a way 
that Px acts the first. As soon as the ordering of operators has been introduced, one 
can work with functions of c-numbers Px and x. These functions are called "symbols" 

mi 



It is shown in Appendix A that leading terms L'^{px^x) of the effective 



Hamiltonians L^{px,x,h) are eigenvalues of H{px,x): 
HiPx,x)xoiPx,x) = L^{px,x)xoiPx,x), 
where Xo{Px,x) are two eigenvectors of the matrix H{px,x] 

Lo{Px,x) = ±\p\+u{x), Xo{Px) = ^ 

We note that in the absence of a magnetic field Xo does not depend on x. The first 
correction Lf{px,x) reads 

dxt dL^ _ ldl4_d(j)p ^ u'{x) Py 
2 dx dpx 




(67) 



(68) 



L^iPx 



X 



dpx dx 



2 pI+pI 



(69) 



Standard semiclassical treatment (see Appendix B ) can be applied to scalar Schrodinger- 
like equations L^ip^ = Eip^ . We are looking for a solution in the form ip^ = 
^is±{x)/hy^±(^^^ /i), A^{x, h) = A^{x) + hAt{x) + ... This gives 



A^ix) 



dpx 



'1/2 



exp 



—i / dx 



, dpx 



-1 



i d'L 



i2 r± 



2 dpxdx 



(70) 




Figure 1. Classical phase space for a) n-p and b) n-p-n junctions. In the electronic 
region the velocity is codirectional with the momentum and in the hole region the 
velocity has an opposite direction to the momentum. Therefore electronic trajectories 
are clockwise oriented but hole trajectories are oriented counterclockwise. Plus and 
minus in figures denote signs of dpx/dx. 



with = dS^ /dx to be found from the Hamihon-Jacobi equation (pi, x) = E, where 



dpx 



\P.\ _ {[E - n{x)]' - Pi) 



1/2 



bl 



\E — u(x) 



Differentiating the Hamilton- Jacobi equation with respect to x we find 



^^0 dpx ^ dL^ 
dpx dx dx 



0, 



whence 



This gives 



dpx 



ip{x) 



1 dL^ 
p'^ dx 

\E-u{x)\^l^ 



1/4' 



^±iS+{x)/h+i,j,^{x)/2 



{E — u{x))^ — pI 
Though it is possible to define locally Xo (EH]) as 



XoiPx] 



-i(f)p/2 



V2 V ±e^<^''/2 



(71) 



(72) 



(73) 



(74) 



(75) 



to obtain Lf = 0, such a choice does not provide a single- valued function in the classical 
phase space. 

Let us first consider a scattering problem for Py ^ following Appendix B and 



Appendix C[ For an electron coming from the left of the classically forbidden region we 

\E-u(x)\'/^ 



have 

■ip{x) 



{E — u{x)Y — pi 



1/4 



^iS+{x)/h+i,j,+ {x)/2 



+ e 



-iS+ {x)/h+i4p {x) /2-i7r/2^ 



(76) 
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and 
where 



\E 



MX, 



11/2 



{E-uix)y-Pl 



1/4 



^i3+{x)/h-i<p+{x)/2 _|_ ^-iS+(x)/h-i(f>-{x)/2-iTT/2 
^iS+(x)/h+icf>+{x)/2 _|_ ^-iS+(x)/h+i(l>p{x)/2-iTr/2 



S 



x] 



± 



xo 



^v'^{x') — p^dx', 
TTSgn{py) 



[X 



Arg (±y^ 



v'^(x) 



iPy 



V[X] 



u[x] 



E. 



(77) 

(78) 
(79) 



Note that 



[x) continuously depends on py when it passes through zero and 0p (x) 
undergoes a jump of 2-n. The reflection coefficient r can be computed from fl76|) or f l77|) . 
It is usually defined as the coefficient in front of the semiclassical solution corresponding 
to the outgoing wave. One can also assume that the potential tends to a constant at 
infinity and take the coefficient in front of the plane wave, which is a particular case of 
the definition given above. Obviously, the refiection coefficient defined in such a way 
does not depend on x. Choosing fl7i|) as incoming and outgoing solutions, we can write 
wavefunctions on the left of the classically forbidden region as 
|E-m(x)|1/2 



ijj^x) 



^ix) 



{E — u{x)y — 
|E-m(x)|i/2 



1/4 



^^iS+ix)/h+i<t>+ix)/2 _^ ^^^^^^-iS+ix)/h+i4>-{x)/2 



0) 



1/4 



n(x))^ 

^iS+{x)/h-i^+{x)/2 



+ r{p, 



y) 



^iS+(x)/h+i4>+(x)/2 

and (CEl), (177]) we conclude that 



g-i5+(x)//i-i(/)p (x)/2 
^~iS+{x)/h+i<t,-{x)/2 



Comparing ( l80l) 

r{py) = e-'^'\ 

A similar calculation for a hole coming from the right gives, see also figure [H 



r{py) 



^m/2 



(82) 



^3) 



We paid attention to the definition of the reflection coefficient, since it may lead to 
discrepancy for the Dirac particle. The problem appears due to a jump of 2tt in 
(j)Z{x) at any fixed X clS 8b function of Py when it goes through zero. This jump is a 
consequence of the cut at 0p = ivr. At any Py this cut corresponds to infinite 
negative x-component of the momentum, and does not imply any discontinuities in the 
region, where the potential is finite. This jump results in the jump of vr in the phase 
of the wavefunction corresponding to the outgoing wave. However, the phase difference 
0p (x)/2 — 0+(x)/2 = TTsgn {py)/2 — (j)p{x) tends to zero when x tends to a turning point 
Xo and can therefore be treated as one half of the angle around the origin in p-space 
accumulating during the motion of a classicle particle from the point x to the turning 
point Xo and back. The peculiar behaviour of the phase difference can mathematicaly 
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be expressed as the noncommutativity of limits: 

lim (a^) - <l>pi^)] = (84) 

The jump in the sign of the outgoing wave must be compensated by a kink in the 
reflection coefficient, since the whole wavefunction should analytically depend on py. To 
get rid of the jump one can redefine the outgoing wave and write ^2] 

\E -U{x)\^/^ riS+(x)/h+i<l,+ (x)/2 _^ ^(^p ^^-iS+(x)/h-i<p+{x)/2\ _ (^g5^) 

\e - u{x)y - Pl]'^' ^ 

Though preserving the analyticity of r{py), such a definition introduces an artificial 
jump of the phase as a function of x upon reflection at negative Py. Therefore we do 
not use it below. 

The reflection coefficient, defined in accordance with (1801) . ( IHTl) does not depend on 
the sign of Py. It is completely defined by the orientation of the phase space, which is 
clockwise for an electron region and counterclockwise for a hole region (see figured]). 
Finally, the reflection can be written as 

r{py) = e^'''l\ (86) 

where '-' corresponds to electron and '+' to hole regions. 

It is important to note, that the phase — 7r/2 and the module 1 of the reflection 
coefficient (!86l) were obtained under the assumption that there is no multiplicity change! 
It is not the case when — )■ and the trajectory in the phase space tends to a separatrix, 
see figure [2j 

Let us now turn to bilayer graphene. The Hamiltonian describing the charge carrier 
dynamics reads 

{px - iVy? 

iPx + iPyY 
Eigenvalues and eigenvectors of H{p^, x) are 



H={ (87) 



We obtain 



L^ = ±p' + u{x), X^ = ^[ ±1 j- (88) 
E — u{x) =Fp^ 

S{x) = r J±[E -u{x)]-pldx. (89) 

Jxo ^ ^ 

Obviously, the result flS^ is valid for the bilayer as well, since the orientation of 
the phase space is the same. 

Between two classically forbidden regions effective Hamiltonians superimpose the 
following quantization conditions (see Appendix C for details): 

^ j p^dx + ^ A0P = 27r + , (90) 
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a) b) 



Figure 2. Classical phase space for n-p-n junction in a) single and h) bilayer 
graphene. Different trajectories correspond to different values of Py. One sees that 
for normal incidence (separatrices) the smoothest classical trajectory corresponds to 
total transmission in single layer graphene and to total reflection in bilayer graphene. 

where /3 = 1, 2 for single and bilayer respectively, i/ = 2 is the Maslov index and A0p 
is the total phase gain along the closed classical trajectory. The term /3A0p/2 is the 
Berry phase in graphene |13]. It is clear that A0p acquires a non-zero value only if 
the trajectory in p-space encloses the origin. Therefore in the absence of magnetic field 
A0p = 0. Quantization condition (!90|) allows one to determine resonance angles. 

Though the considered diagonalization is very powerful to deal with complicated 
matrix Hamiltonians in a classically allowed region, it possesses a substantial 
disadvantage: it treats electrons and holes separately neglecting tunneling effects. In the 
classically forbidden region when \p\ = 0, i.e. px = iPy effective Hamiltonians Lq become 
degenerate. At this point electron to hole transition may occur and the diagonalization 
fails. This transition is the origin of the Klein tunneling. 



8. Normal incidence 



In the case of normal incidence Py = and at the point xq, where u{xo) = E there is a 
multiplicity change, i.e. effective Hamiltonians Lq become degenerate (see figure [2]). To 
study wavefunctions in this case one can not apply a standard semiclassical treatment, 
described in Section [71 since there may be a "jump" between L"*" and L~ . Fortunately, 
for the normal incidence in graphene there is an exact pseudospin conservation, which 
allows one to study this case in detail. 
For = equations ( l62l) . ( j63l) read 

[a^pf + m(x)-E]^ = 0, (91) 

where /3 = 1, 2 for single and bilayer respectively. Eigenvectors of a^p^ do not depend 
on , therefore ( l9T|) can easily be diagonalized, which leads to 

[±pf + u{x) - E]r],,2 = 0, (92) 

where 

^ = ( J ) ^1 + ( _\ ) ^2. (93) 
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Figure 3. Comparison of tlie initial potential and the real part of the effective 
potential. 

In this case the eigenvalue of ax ( "pseudospin" ) persists. Pseudospin conservation leads 
to very different physical consequences for single and bilayer graphene. 

For single layer graphene pseudospin conservation means the conservation of the 
x-component of the velocity. Equation f l92|) is the first order differential equation, which 
can be solved exacly. We obtain 

771,2 = Ci,2 exp [±i £ [E - u{x')]dx'^ , (94) 

where Ci,2 are some constants. The absence of the reflected wave in flM|) means that for 
any potential shape one has a perfect transmission. Thus we conclude that at the point 
Px = there is a total transition between electron and hole states since Hamiltonians 
f l68|) depend on \p\ in contrast to (!92|) ! 

For bilayer graphene pseudospin conservation is equivalent to the conservation of 
particle type, as is seen from the comparison of flHHj) and f l92|) . Therefore, an incoming 
particle obeys the Schrodinger equation (1921) everywhere. For a "Klein-setup" this leads 
to exponentially decaying transmission as a function of a potential width and height. 

Total transmission for normally incident electrons in single layer graphene and its 
exponential damped behaviour in bilayer have natural explanations of classical phase 
space (figure |2]). In both cases the most probable process corresponds to the smoothest 
trajectory, constructed from separatrix pieces. For the single layer such a trajectory goes 
through the barrier and gives total transmission, while for bilayer one has to choose the 
trajectory reflected from the barrier to avoid discontinuity in the second derivative. 

9. Exact reduction to effective Schrodinger equations 

In Section [7] it was discussed that the standard adiabatic diagonalization fails to describe 
Klein tunneling, since it treats electrons and holes separately. However the existence of 
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the exact diagonalization (Section [H]) for a normal incidence raises the issue of a possible 
generalization for angular scattering. We shall show that for a single layer graphene 
there exists an exact transformation, reducing the original Dirac equation to a scalar 
Schrodinger-like equation with a complex potential. It is clear that such a procedure 
can not be a unitary transformation of the original Hermitian Hamiltonian. 

Let us turn back to the Dirac equation for single layer and write it in the form 

(o-p + t;(x))^ = 0, (95) 

where p = {px,Py), ^{x) = u{x) — E. Let us act on the last equation from the left by 
the operator crp — v{x). Then we get 

(crp - v{x)){ap + v{x))'^ = {pl+pl - v{xf + a^[p^,v{x)])-<^ 

= iPl +pI- v{xf - iha^v\x))'^ = 0. (96) 

Remarkably, ( l96l) contains only the single matrix ax- Therefore it can easily be 
diagonalized. We write 

^=(]]vi+( \]v2 (97) 



and obtain 



1 / \ -1 



h^^ + v{x)^ -pl±thv'{x)]r]i^2 = 0. (9^ 



Functions 771^2 are not independent and the connection formula can be obtained from 
(1951) . Function 772 can be reconstructed from rji as 

T]2 = — ( h-^ + iv{x)] rji. (99) 
Py \ dx J 

In figure [3] one sees the initial potential landscape and the real part of the effective 
potential in f l98|) . Though this equation takes the form of a Schrodinger equation, there 
are two substantial distinctions as compared to a common Schrodinger particle: i) the 
effective potential is complex and ii) it depends on energy. 

10. Single n-p junction 

10.1. Exact solution in the case of a linear n-p junction 

Let us first consider an exacly solvable model for a linear potential v{x) = ax, a > 
|44] . Introducing a new variable x' = (a/hY^'^x and new y-component of the momentum 
Py = {ha)~^^'^Py we exclude h and a from ( 198|) . Then it takes the form (we omit primes): 

(^+x^-pj + z)r/i=0. (100) 



Introducing a new variable z such that x = we have: 
V dz 



2 + e'(^-p') + evUi = o. (101) 
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Choosing an appropriate value for ^ we can reduce fllUlj] to the Weber's equation 



(102) 



Indeed, choosing ^ = e ^"^^^f \/2 and solving the Weber's equation |l5lll6] we obtain 

r|^ = c^D.i^Ple'^l^x) + C2/^-.-i(V2e3'"/^x), (103) 

where v = ipy/2 and D,y are the parabolic cylinder functions. For these functions the 
following identities hold: 



= ^D^_i{z) - -D^{z), = i^D^[z) - D^+i[z) 



t/4 



Applying the first equality from (11041) we find 

+ i:^ D,{y/2e'^'S) = ^ve'^I^D,^^{^e'^ 
From the second equality in (11041) we have 

[^^^^ D_,_i(v^e=^^'^/^x) = V2e-^"/^D_,(y2e=^*"/^ 
Substituting ffTOSjl . ffT06|) into ([99]) we obtain: 



(104) 



(105) 



(106) 



^2 



Py \dX 

Cl 



+ ix\ r]i 



y2z/e*"/^D^_i(y2e^"/S) + —V2e-'''/^D_^{y/2e^''' 

Py Py 



/4 



X) 



(107) 



From f lT03|l . ( HMf we have for ^ = (^i,^/'2): 

^1,2(2;) = mix) ±^72(3;) 

= Cl ( D,{V2e'^/^x) ± 



Py 



D,^i{V2e'^/^x) 



+ C2 D-.-i(v^e=^^-/^x) ± ^ D_,(v^e3-/^x) 



108) 



Using the asymptotic expansions of the parabolic cylinder functions (see 



Appendix D ), we find when x — > 00 (hole region): 



+ C2 



^27r 



V -ix^ l2-i-K{v+\) _|_ 

r(z/ + i) 2 



1) ^ V2e-*-/4 



2:0 e 



Py 



(109) 



and when a; — )■ — 00 (electron region): 



V'1,2 Cl 



^2 



r(i - U) Py 



± (%)-V^ 

Py 

where zi = \/2e^'^^^\x\, Z2 = \/2e^^'^^^\x\ and a bar means complex conjugation. 



;iio) 
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Now we turn to the discussion of the scattering problem. While tunneling through 
the barrier the Dirac particle turns from an electron to a hole or vice versa. The 
x-component of the group velocity of the hole Vx = OLq /dpx = —Px/p has an opposite 
sign with respect to its momentum p^. Let us consider an electron, coming from 
—oo with a positive velocity Vx- It corresponds to the action S^{x) ~ —x^/2, since 
Px = dS'^{x)/dx ~ —X > and Vx = Px/p > 0. Thus, the reflected electron 
corresponds to S^{x) ~ x'^/2. The transmitted hole with a positive velocity has a 
negative momentum px- Hence it corresponds to the action S^{x) ^ — x^/2. From the 
absence of the incoming wave in the hole region we find C2 = 0. 

Let us consider (IHTl) at infinity. Then we have for the action: 



S+(x) 



sgn {x)\py\ 



-sgn [X 



1/2 - pMy 




pi In 



X 

Py 



\ 




where we assumed that |x| > \py\. For large x we obtain 



S+{x) 

Thus for large negative x (ISTj) reads 

= ^-ix^/2+ipl/i+{i/2)pl\n(2\x\/\py\) 



^sgn (x) - ^ - pjln ^ 




111) 



;il2) 




+r(pj^)e*'''/^"*^^/'^"*^*/^^^« in(2|x|/|py|) 



g-i7rsgn(p„)/2 
giTT sgn(py)/2 



This gives for the reflection 



2 /o tJ2 /o^ ^^IJ2 



(113) 



(114) 



where 6{py) = Py/2 — {Py/2) \n{py/2) — tx /A. Using equalities 



r(i-z/)r(i 



Sm TTZ/ 



ttpI 



we can write the reflection coefficient as 
r{py) = \fr- 



(115) 



;ii6) 



where 7(py) = Argr(l — ipy/2). From the asymptotic expansion of the F-function at 
large arguments [15], one concludes that 7(py) tends to 0{py) when py tends to infinity. 
At small Py the reflection coefficient is proportional to \py\. This nonanalytic behaviour is 
compensated by the jump of the phase of the reflected wave as we discussed in Section [3 
Comparing the coefficient in front of incoming and transmitted waves we find for 
the transmission amplitude 

(117) 



t 



gi7rsgn(py)/2gi7ri^ 



gi7rsgn(pa)/2g-7rp2/2^ 
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For the transmission probability we thus have 



-7rp^ sin^ <f>p 



:ii8) 



This result was first obtained by Cheianov and Fal'ko [S]. Considering the scattering 
from the right to the left we find that the transmission amplitude in this case is 



(119) 



The transfer matrix connecting incoming and outgoing waves from the right to the 
left of the barrier for positive a is 



rp _ ^-iTTSgn{py)/2 



1/2 



(g-p? - ly 



,i{7-e-7r/2) 



For negative a the transfer matrix reads 
/ 



J'_ — g*'rsgn{py)/2 



(fp'^y - l) 



1/2 



(f^'^y - l) 



1/2 



-7-71/2) 



,i(7-9-7r/2) 



3 7rp2/2 



(120) 



;i2i) 



10.2. Transmission probability in semiclasical approximation 

Let us consider an outgoing hole on the right of a generic potential monotonously 
growing from the left to the right. In semiclassical approximation it is described by 
the wavefunction 

g-*.sgn(p,)/2|^_^^^^|l/2 / g*<A?W/2 

' g-i</'+{:r)/2 



{E-u{x)y-pl 



1/4 



(122) 



Using equalities 



cos 



COS I 



2 



1 + cos < 



' 1 — COS ( 



sm 



\Px 



we write it as 



-in sgn (py)/2 



2 ' V 2 ^ 
[{E ~ ujx))' - pI]'/' 
\E - u(x)\ 



iS+/h 



1/4 



V + . V 



Pi 



1/2' 




(123) 



(124) 



with v{x) = u{x) — E. According to f ll24p the components of '^{x) can be represented 
exactly as a sum or a difference of functions rji, rj2 obeying Schrodinger-like equations 
( l98l) with a complex potential. Despite the complexity of the potential, according to 
to connect a transmitted wave on the right of the barrier with an incoming wave on the 
left of the barrier one can still use an analytic continuation in the classically forbidden 
region similar to |18]. In contrast to a usual Schrodinger equation the transmitted 
hole has a negative momentum, so \E'(x) in f ll22p is proportional to g-^'S'^/'^^ not to 
giS+//i_ Therefore the passage should be done in the lower complex half-plane. Since 
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both functions t]i and 772 allow analytic continuations in the lower half-plane, can 
also be continued in the lower half-plane. 

Performing the passage and connecting the outgoing wave with an incoming wave, 
we obtain for the transmission coefficient: 

fX2 



-j- _ ^iiTSgn{py)/2^-K/h 



Xl 



\Jp1 ^ v'^{x)dx 



(125) 



where Xi and X2 are two turning points, i.e. solutions for the equation Py — v {x) = 0. 
The standard complex WKB technique |l9l |50l ST] does not allow one to compute 
the reflection coefficient with an exponential accuracy. The module of the reflection 
coefficient can be computed from the unitarity of the scattering matrix: |rp + |tp = l. 
This gives 



l-e-2^/^ (126) 

The semiclassical phase of the reflection coefficient can be reconstructed from 
Thus in semiclassical approxiamtion we obtain: 



- e-2^/'^e^'"/^ (127) 

where '-' corresponds to the electron and '+' to the hole region. For small Py any 
potential can be linearized in the classically forbidden region. Comparing flllTp and 
(11251) we see that in the limit Py ^ the semiclassical transmission becomes exact. 
Therefore (11251) can be used as a uniform approximation for the transmission coefficient 
at any py. Then the uniform approximation for the reflection coefficient reads 



1 _ e~^K/h^Tirr/2+ie^ (^28) 

where tends to zero when py tends to infinity. For the phase we used the expression 
[511 152] 

^ K K ^ f K\ TT iK\ 
which was obtained by the replacement vrp^/2 by K/h in 6{py) — l{Py)- 



11. Klein tunneling in n-p-n junctions. Fabry-Perot interferometer 

Let us consider some generic potential barrier u{x). In graphene it implies n-p-n junction 
(see figure [3]). In terms of fl98l) n-p-n junction becomes a complex double-hump potential. 
The transfer matrix in this case is 

/ iS/h Q \ .X3 I 

T = T^[^ ^-^s/h)T-, S = j^^ ^v\x)-pldx, (130) 

where we assume that X2 < X3. There is no extra phase coming from 0^ since each of 
these functions takes the same values at both turning points X2, Xs, lying in the hole 
region. The transmission coefficient reads 

1 p-iS/hp-K^/h-K2/h 
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The obtained transmission amplitude can easily be treated in terms of a sum of 
probability amplitudes of multiscattering processes leading to transmission |12]. One 
sees that for normal incidence Ki = K2 = and the module of the transmission 
coefficient becomes one. Transmission resonances can be found from the condition: 

which coincides with the quantization condition (!90|) for large py. For a symmetric n-p-n 
junction the resonant transmission is always one, since Ki = K2. For an asymmetric 
junction resonant values of transmission decay as 

1 

^ ^ cosh{Ki/h - K2/h) ^^^^^ 
when Ki/h ^ 1 and K2/h ^ 1. From (11331) one sees that resonant values of 
transmission exponentially decay as a function of \Ki — K2\. Such a fast decay can 
be crucial if one wants to weaken the influence of side resonances. 



12. Numerical results 



In figures HUE] we compare our semiclassical predictions with numerical results, obtained 
from a multistep approximation of the initial potential. A check on the accuracy of the 
calculation for constant Py is provided by the current 

J, = ^V,^, djjdx = 0. (134) 

To simulate an n-p junction we used the potential 

V{x/h) = 0.5 U^ax [1 + tanh(10x//i - 5)] (135) 

with a characteristic length scale li. An n-p-n junction was simulated as an n-p junction 
with a characteristic length /i, a p-n junction with a characteristic length I3 and a 
constant potential in between of the length I2. 

In figure m we show the comparison of our numerical result for an n-p junction with 
the semiclassical transmission (I125P and the transmission for a linear potential (I117p . 
While the semiclassical prediction works uniformly over the entire range of angles, the 
prediction obtained from a linear potential works only for small angles. 

In figure[5]we show the comparison of our numerical results with the prediction (11310 
for a symmetric n-p-n junction. We also show the semiclassical result, which is obtained 
by setting Bi = 62 = 0. The agreement between the latter answer and numerics gets 
better as the angle increases, i.e. deep in the semiclassical regime. The result (11251) 
uniformly approximates the numerical data over the entire range of angles. 

In figure [6] the result for an asymmetric n-p-n junction is shown. The height of 
resonances is seen to decay. The suppression of side resonances for asymmetric junctions 
in single layer graphene can have essential consequences for attemps to confine Dirac 
particles! 

Numerical computations for bilayer graphene using the above procedure are less 
accurate, due to the presence of real exponentials everywhere. Therefore we were unable 
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Figure 4. The angular dependence of the transmission coefficient for a particle of 
energy 80 meV incident on an n-p junction of height 200 meV. The potential is given 
by equation (|135|) . with li = 70 nm. The blue line shows the numerical result with 49 
steps, the dashed line shows the semiclassical result ()125|) and the red line shows the 
result for a linear potential ()117|) . where the parameter a was taken as the derivative 
at the central point of the junction. 




Figure 5. The angular dependence of the transmission coefficient for a particle of 
energy 80 meV incident on an n-p-n junction of height 200 meV. The barrier width 
I2 — 250 nm and n-p and p-n regions have characteristic lengths h = I3 = 100 nm. 
The blue line shows the numerical result for 99 steps, the red line shows the uniform 
approximation (|13ip and the orange line shows the semiclassical answer {Qi = 82 = 0). 
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to check the quantization condition fl90p numerically with a high precision. To check 
the accuracy of the computation we used the current 



In figure [7] we show our numerical results for a symmetric and an asymmetric n-p-n 
junction with the same shape of the potential as before. In contrast to the case of single 
layer, resonances do not seem to decay in this case. 

13. Conclusion 

Let us summarize our main results. The detailed analysis of the refiection-and- 
transmission problem for the Dirac electrons demonstrates essential differences from 
the conventional Schrodinger case, due to the role of the Berry phase. The reflection 
coefficient turns out to be nonanalytic function of the transverse momentum Py vanishing 
as \py\ at Py —J- 0. 

We have presented a complete treatment of the chiral tunneling for both single and 
bilayer graphene in terms of a classical phase space. This gives a natural explanation of 
complete transmission of normally incident wave for single layer and its exponentially 
damped transmission in bilayer. We have also demonstrated that, for the case of 
nonsymmetric n-p-n junction in single layer graphene, there is total transmission for 
the normal incidence only, and other maxima are suppressed. Our numerical studies 
show that for the case of bilayer there are always magic angles with total transmission. 
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Appendix A. Effective Hamiltonians in adiabatic approximation 

In this part of the Appendix we show how to reduce in adiabatic approximation an 
initial matrix Hamiltonian to a set of effective scalar Hamiltonians. Our consideration 



Let us consider an eigenvalue problem for an Hermitian matrix Hamiltonian H, 




(136) 



(137) 



follows [all M 





Chiral tunneling in single and bilayer graphene 27 




Figure 6. The angular dependence of the transmission coefficient for a particle of 
energy 80 meV incident on an n-p-n junction of height 200 meV. The barrier width 

12 = 250 nm and the n-p and p-n regions have characteristic lengths h = 150 nm and 

13 = 50 nm, respectively. The blue line shows the numerical results for 99 steps, while 
the red line shows the uniform approximation (|131|) . 




Figure 7. The angular dependence of the transmission coefficent for a particle of 
energy 17 meV incident on symmetric and asymmetric n-p-n junctions in bilayer 
graphene. Each junction has a height of 50 meV and a width I2 = 100 nm. The 
blue line shows the numerical result for a symmetric junction with h = I3 = 10 nm, 
while the red line shows an asymmetric junction with li = 20 nm and I3 = 40 nm. All 
calculations were done with 99 steps per junction. 
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where we assume that = —ihd/dx acts first and x acts second. In what follows we 
always assume this operator ordering. Let us introduce a vector operator x and a scalar 
wave function ip by the requirement 

^{x) = xi-^hd/dx,x,h)ij{x) . (A. 2) 

We want to satisfy an eigenvalue problem 

L{~ihd/dx,x,h)ip{x) = Eip{x) , (A. 3) 

with an effective Hermitian Hamiltonian L. Substituting flA.SP in flA.ll) we obtain: 

m-xLMx) = 0. (A.4) 
The last equality will be fulfilled for any iIj{x) if the following operator equality holds: 
H{— ihd/dx, x)x{— ihd/dx, x,h) = x{— ihd/dx, x,h)L{— ihd/dx, x,h). (A. 5) 

We will solve it by passing to symbols of operators (see [iQl 1^): 

smb[A{— ihd/dx, x)B{— ihd/dx, x)] = A{px — ihd/dx, x)B{px,x) . (A. 6) 
Applying this formula to the above case, we obtain 

H{pj. — ihd/dx, x)x{Px,x,h) = x{Px — ihd/dx, x,h)L{px,x,h) . (A. 7) 

Let us expand this expression with respect to the parameter h ^ 1. In zeroth order we 
get 

HiPx, x)xoiPx, x) = Lo{p^, x)xoiPx, x) , (A.8) 

where Lo{px,x) = L{px,x,0), Xo{Px,x) = x{Px,x,0). The first order term in the 
expansion gives 

~ ^ + = ^ + ^oXi + XoLi , (A.9) 

OPx OX OPx OX 

where Li and xi ^i-re the first order terms in L and x with respect to h. The above 
expression can be rewritten as 

[H-Lq)xi=i^ T. — + xoLi. (A. 10) 

OPx OX OPx OX 

Let us multiply the last equation by Xo from the left. Since both H and Lq are Hermitian 
matrices, we obtain 

dx ^ dpx dx 
where we used the equality XoXo = 1- 
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Appendix B. Semiclassical approximation 

Appendix B.l. x -representation 

To solve equation flA.SP we will use the semiclassical ansatz 

^{x)=e'^^''^'^A{x,h) . (B.l) 
Then (]A.3[) can then be rewritten as 

L{dS/dx - ihd/dx, x, h)A{x, h) = EA{x, h) . (B.2) 

This equation can be expanded order by order in x which to zeroth order gives the 
Hamilton- Jacobi equation 

Lo{dS/dx,x) = E . (B.3) 

From this equation we can determine the action S{x), as is well known from classical 
mechanics [53]. To first order in h we obtain the equation 

i d^Lo d^S 



.dLodAo 
- I- h LiAq 



;Ao = Q 



(B.4) 



dpx dx 2 dpi dx"^ 

where all terms should be evaluated at p^ = dS/dx. When multiplying by the amplitude 
Ao, equation (]B.4p can be rewritten as 



i^d_ ( ^Lo^2 
2dx\ dpx ° 



A =0 



(B.5) 



2 dpxdx ^ 

where the total derivative acts on both x and px = dS/dx. This equation can be solved 
exactly to determine the amplitude 

-1 



An 







-1/2 




dpx 


exp 




dLo 


-1/2 




dpx 


exp 



dx 



dx 
\dpx 



dpx 



2 dpxdx 



M 



where we have defined M. Using equation ( 1A.11I) it can be written as 
M 



^dHdxo , ^dxodLo 1 O'^Lq 



dpx dx dpx dx 2 dpxdx 
Appendix B.2. p -representation 

We can also solve equation (IA.3P by passing to p- representation [iQldl]: 



(B.6) 
(B.7) 

(B.8) 



(B.9) 



L{px,ihd/dpx,h)ilj{px) = E^{px) . 

In this equation p^ acts first, while ihd/dpx acts second, contrary to the case we 
considered before. To solve this equation, we use the semiclassical ansatz 

^{Px) = e-''^^-^/''A{p,,h). (B.IO) 
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Similarly to equation (IB.2p equation (IB.9P can be rewritten as 
L{p^, dS/dp^ + ihd/dp^, h)A{p^, h) = EA{p^, h) , 
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(B.ll) 



When this equation is expanded order to order in /i, one obtains to zeroth order the 
Hamilton- Jacobi equation 

Loip.,dS/dp,) = E , (B.12) 

from which the action S{px) can be determined. The first order term becomes 



^— — ^ -^1^0 + o^rr^-TT^o + ^TrrTT— = , 



(B.13) 



dx dpx 2 dx'^ dpi dxdpx 

where all terms have to be evaluated at x = dS/dp^- After multiplication by Aq one 
finds 



id/ dL 
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2 dpx \ dx 



i d^Ln 



2 dpxdx 



which can be solved exactly to give 



An 





dLo 


-1/2 




dx 


exp 




dLo 


-1/2 




dx 


exp 



dpx 



dpx 



'dLo' 
'dLo 



-1 



i^i + 



1 d^Lo ' 

2 dpxdx 



-1 



dx 



M 



Using equation flA.llI) the quantity M can then be written as 

M = -xo — — + Xo — ^ + . 
dpx dx dpx dx 2 dp^dx 



(B.14) 

(B.15) 
(B.16) 

(B.17) 



Appendix B.3. Matching 

Since the solutions (IB.ip and fIB.lOp come from the same equation, they should be 
related. To find out what this relation is, we look at the Fourier representation of ( IB. II) . 
which is defined by 

V2nh 

Since the parameter h is assumed to be small, we can calculate this integral using the 
stationary phase method. The result is 



(B.18) 



Ao{xs) 



J{S{xs)-PxXs)/h+i sgn(S"(xs))7r/4 



(B.19) 



\S"{Xs)\ 

where the point Xg = Xs{px) at which the phase is stationary is to be found from the 
equality 

S'{xs)=Px. (B.20) 
From the comparison of fIB.lOp and (IB.19P we find 



^isgn{S"{Xs))7T 



Px 



(B.21) 
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Figure CI. The phase space of a classical particle is covered by maps /, //, ///, 
IV. Regular maps / and III can be uniquely projected onto the x-axis, while singular 
maps // and IV can be uniquely projected onto the px-Bxis. The maps are chosen to 
overlap. 



Appendix C. Bohr-Sommerfeld quantization rule 

Now let us consider the phase space which is shown in figure ICll The part / of the 
phase trajectory can be projected onto the x-axis. Therefore we can use a standard 
WKB ansatz in this region. On the contrary, in region II we can use a standard 
WKB ansatz in p-representation. Since regions / and II overlap, functions are related 
according to ( 1B.21[) . Since S"{x2) < 0, we have ipi{x) — )■ 

It is easily seen that the above reasoning can also be applied to regions II and ///. 
Since S"{x3) > we obtain ipiii{x) — )• '^'nipx)^™^'^ , or 

iJi{x)^ijju{x)e-'^/^ . (C.l) 

So, passing the turning point lying in the region // we have picked up an extra factor 
exp(— ?7r/2). Since S"'(x4) < and S"{xi) > we pick up another factor of exp(— ^72) 
when we go through region IV and pass the second turning point, 

i^iiAx) ^ Mx)e-'^'^ . (C.2) 

In passing one full turn along the circle, one therefore sees that ipilx) — il)i{x)e~'^'^ . 
The wavefunction should be single-valued, which means that the exponent should be a 
multiple of 27r. From equations ( IB. II) and ( IB. 71 ) we therefore find 

!)B - TT = 27rn , (C.3) 



where px = S'{x) is to be found from the Hamilton- Jacobi equation ( IB. 31) . The quantity 
is the Berry phase, defined by 




^dII_dxo^ t dxo dip ^ 1 d'^Lp \ 
dpx dx dpx dx 2 dpdx J 
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Equation (IC.Sp can be rewritten as the Bohr- Sommerf eld quantization rule 

Looking back at the above derivation one sees that the term 1/2 can be written as z//4, 
where u is the number of turning points (Maslov index in this particular case) [H]. 

Appendix D. Asymptotic expansions of parabolic cylinder functions in 
different Stokes sectors 

For completeness we placed in this section the asymptotic expansions of the parabolic 
cylinder functions in different Stokes sectors at \z\ — t- oo according to 

r e-^'/S^(l + 0[z-^]), -7r/2 < arg (z) < 7r/2, 

DAz)^\ e-^/S-(l + 0[z-^]) - e-^/^--v^^— - (1 + o[z-^]), arg(z)<-7r/2 (D.l) 

[ e-^'/S^(l + 0[Z-']) - e£^li^^rvj^(i ^ ^^^-2])^ > 

We assume that — vr < arg (2;) < n. 
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